
# This file creates Figure SI.2 -------------------------------------------

# install.packages("rmeta")
# install.packages("ggplot2")
# install.packages("dplyr")

library(rmeta)
library(ggplot2)
library(dplyr)

#read the estimates data
robust_coefs <- 
  read.table("robust_coefs.txt")

# assign variable names
vars <- unique(robust_coefs$labe)

# create empty data frame
meta_effects <- 
  data_frame()

# Find and assign meta estimates
for(i in 1:length(vars)){
  for(j in 1:2){
    FE <- ifelse(j == 1, "YES", "NO")
    var_i <- vars[i]
    
    var_data_mother <- 
      robust_coefs %>% 
      filter(labe == var_i & parent == "Mother" & fe == FE) 
    
    var_data_father <- 
      robust_coefs %>% 
      filter(labe == var_i & parent == "Father" & fe == FE)
    
    meta_mother <- 
      meta.summaries(var_data_mother$coef, 
                     var_data_mother$se)
    
    meta_father <- 
      meta.summaries(var_data_father$coef, 
                     var_data_father$se)
    
    meta_effects <-
      bind_rows(meta_effects, 
                data_frame(var    = var_i,
                           effect = c(meta_mother$summary,
                                      meta_father$summary),
                           se     = c(meta_mother$se.summary,
                                      meta_father$se.summary),
                           parent = c("Mother", "Father"),
                           place  = i + (j - 1.5)/5,
                           fe     = FE))
  }
}

# plot meta estimates

plot_meta <- 
  ggplot(data = meta_effects, 
         aes(x = effect, 
             y = place,
             xmin = effect + qnorm(0.025) * se,
             xmax = effect + qnorm(0.975) * se,
             alpha = fe)) + 
  geom_point(size = 4) +
  facet_grid( ~ parent) + 
  geom_errorbarh(height = 0, size = 2) + 
  scale_y_continuous("", breaks = unique(round(meta_effects$place)),
                     labels = vars) + 
  theme_classic() + 
  scale_x_continuous("Twinning vs non-twinning") + 
  geom_vline(xintercept = 0, alpha = 0.5, linetype = 2) + 
  scale_alpha_discrete("Age FEs", range = c(0.5, 1)) + 
  theme(legend.position="top") + 
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 16),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        title = element_text(size = 16),
        panel.spacing = unit(1, "lines")) 
